#!perl

use strict;
use Getopt::Long;
use MaterialsScript qw(:all);


my $doc= $Documents{"Mg-water.xtd"};

my $file1 = Documents->New("Hbonds_time_vs_number.txt");
my $file2 = Documents->New("Hbonds_distance_vs_number.txt");
my $file3 = Documents->New("Hbonds_O-H-O_distance_vs_number.txt");

my $start = 12;
my $end = 40;
my $step = 0.2;

#------------------------------------------------
my $NumberOfMolecules = $doc->UnitCell->Molecules->Count;
my $trajectory = $doc->Trajectory;

for (my $frame=1; $frame<=$trajectory->NumFrames; ++$frame){
	$trajectory->CurrentFrame = $frame;
	my $newfile = Documents->New("temp_file.xsd");
	$newfile->CopyFrom($doc);
	$newfile->CalculateHBonds;
	my $HBonds = $newfile->UnitCell->HydrogenBonds;
	
	my $NumberofHBonds = $HBonds->Count;
	$file1->Append(sprintf "%d,%d\n", $frame,$NumberofHBonds);	


	$newfile->Delete;
}

	
for (my $i = $start; $i < $end; $i = $i+$step){
	my $j = 0;
	my $k = 0;
	for (my $frame=1; $frame<=$trajectory->NumFrames; ++$frame){
		$trajectory->CurrentFrame = $frame;
		my $newfile = Documents->New("temp_file.xsd");
		$newfile->CopyFrom($doc);
		$newfile->CalculateHBonds;
		my $HBonds = $newfile->UnitCell->HydrogenBonds;
		foreach my $HBond(@$HBonds){
			my $z = $HBond->MidPoint->Z;
			if (( $z > $i) and ($z < ($i+$step))){
				$j = $j +1;
			}
			if ( $HBond->Acceptor->ElementSymbol ne "Cl"){
				if (($z > $i) and ($z < ($i + $step))){
					$k = $k + 1;
				}
			}
		}

		$newfile->Delete;
	}
	$file2->Append(sprintf "%f A,%d\n", $i, $j);
	$file3->Append(sprintf "%f A,%d\n", $i, $k);

}
	

$file1->Close;
$file2->Close;
$file3->Close;


